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ABSTRACT 

This  report  describes  a  numerical  method  to  calculate  the  flow 

of  a  compressible  fluid  in  the  presence  of  shocks.  It  is  related  to 

the  method  of  von  Neumann  and  Richtmyer^^  and,  to  a  lesser  extent,  to 

(2) 

that  of  P.  Lax.  This  method  permits  a  number  of  variations  which 
were  compared  by  testing  them  on  some  analytically  known  one-dimensional 
flow  patterns.  The  calculations  were  carried  out  on  the  Los  Alamos 
"MANIAC." 
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1 .  General  Remarks 


The  one -dimensional  hydrodynamic  equations  can  be  written 
in  the  Lagrangian  form 


and 


where 


Po  -AH  = 

r  dt 

-  , 

(1) 

-H  - 

>  u 

(2) 

^  p  _ 
bt 

-  <  Pow)2  ir  - 

-p  wa  i 

to  a* 

(3) 

u  =  velocity, 
p  =  pressure, 
v  =  specific  volume, 

Po  =  initial  density, 
x  =  initial  position  of  a  particle, 

and  where 

in  which 


S  =  entropy. 

The  first  two  of  these  equations  are  in  the  form  of  conservation 
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theorems  and  can,  with  the  aid  of  Green's  theorem,  be  integrated  in 
a  shock  region.  The  third  equation  is  clearly  wrong  in  the  presence 
of  a  shock.  Nevertheless,  we  use  it  to  expand  p  for  use  in  eq. (l): 

P  -  p(x,o)  -pov2  -hi  t  .  (4) 

Upon  integrating  eq.(l)  to  a  small  value  of  t  =  &  t,  we  get 

u(x,  £t)  -  u(x,0)  -  [  p(*>0)  -  ^§—  U  ]  St  •  (5) 

An  alternative  method  of  integration  is  due  to  Riemann,  who  defined  a 
quantity 


and  showed  that  the  quantities  u  +  CT  remain  constant  along  the 
so-called  characteristics  defined  as  the  lines 


dx 

dt 


=  ±  w. 


(7) 


Therefore,  if  xx  and  xg>  xl  are  two  points  on  the  x-axis  and  if  the 
forward  characteristics  through  x^  and  the  backward  characteristics 
through  x2  meet  a  little  time  later  at  a  point  x,  we  have  at  that 
point 
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(8) 


u  +  cr  =  +  <r 

and 

u  -  O'  =  u2  -  a- 2  ,  (9) 

which  yield 

u  =  |  (  ux  +  u2  )  +  |  (  tr1  -  a~2  ).  (10) 

This  solution  is  related  to  the  construction  of  solutions  of  the  wave 
equation  by  Huygen’s  principle.  This  principle*  states  that  a 
solution  in  a  space  of  an  odd  number  of  dimensions  depends  only  on 
the  starting  values  at  the  border  of  the  domain  of  dependence.  In 
an  even  dimensional  space,  it  depends  also  on  the  interior  of  this 
domain  although  the  influence  of  the  border  is  still  the  stronger. 

In  no  case  do  points  outside  of  the  domain  of  dependence  have  any 
effect. 


*  A  detailed  discussion  can  be  found  in  (3),  Chapter  6,  para.  5.3. 
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2 .  Development  of  Method 


The  computing  scheme  proposed  here  uses  a  rectangular  lat¬ 


tice  (x  ,t  ).  Velocities  u.n  and  displacements  X. “  are  defined  at 
'in  i  x 


the  lattice  points,  whereas  specific  volumes  v ^  +  and  Pressures 


Pi  +  1/2 


are  defined  between  them  as  indicated  by  the  i  +  l/2. 


n  1  , 

If  eq.(lO)  is  applied  to  calculate  u  ,  the  points  1 


and  2  lie  as  indicated  in  Fig.  1  and' can  be  computed  from  eq.(7) 


Fig.  1 

One  can  show  that  the  second  term  in  eq.(lO),  if  expressed  in  terms 
of  quantities  defined  in  the  latter  scheme,  is 


°~i  igk  g  (^i  -  i/g :  pi  ±  l/g*  $t  mi 

2  '  fo^i  - 1  -  X1  -  1J 

This  expression  is  the  same  as  one  would  get  from  the  pressure 
derivative  term  in  eq.(5)-  The  first  term  in  eq.(lO)  can  be  ex¬ 
pressed  in  more  than  one  way.  The  choice  is  introduced  by  using 
different  interpolation  schemes  for  expressing  u1  and  u£  .  If  one 
fits  u  at  the  three  points  x^  _  x^,  and  x^  +  by  a  parabola 

which  has  its  axis  parallel  to  the  u-axis  and  if  one  fits  w  by  a 
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where 


n  1 

qi  +  1/2  "  2 


n 

+  1/2 


/  n  n  v 
(u.  -  u.  -  ) 

i  i  +  1' 


X.  -  -  X. 

1+1  1 


8  t  . 


(12) 

(13) 


This  is  the  same  expression  one  would  get  from  eq.(5)>  If  one  fits 
u  by  two  straight  lines  between  adjacent  points  and  w  in  each 
interval  by  its  central  value ,  one  obtains  eq.(l2)  but  with  a 
different  q: 


n  1  n  /  \ 

qi  +  1/2  "  2  Co  Wi  +  1/2  W  "  Ui  +  1^  ' 


One  might  think  of  other  combinations  of  interpolating  the  two 
quantities  u  and  w,  but  the  two  which  were  chosen  have  the  virtue  of 
leading  to  an  equation  of  the  type  (12),  which  ensures  the  conserva¬ 
tion  of  momentum. 

The  transcription  of  eq.  (2),  which  expresses  the  conserva¬ 
tion  of  mass,  can  be  done  as  follows.  We  introduce  the  position 
X^n  and  integrate: 


.  n  +  1 
i 


n  1  /  n  +  1 

l  2  N  i 


(15) 
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and  calculate  the  volume  from 


n 

i  +  1 


(16) 


As  was  pointed  out  earlier,  eq.(3)  is  in  the  form  of  a 
conservation  theorem.  In  order  to  put  the  conservation  of  energy 
into  evidence,  we  write 


Po  £  <E  *  r>  - 


(IT) 


where  E(p,v)  is  the  internal  energy.  One  has  to  exercise  care  in 
performing  averages  in  the  transcription  of  this  equation  into 
difference  equation  form,  for  u,  which  occurs  on  both  sides  of  eq.  (17) 
is  not  defined  at  the  half  points  as  are  p  and  E.  It  turns  out  that 
one  can  obtain  a  much  simpler  form  if  one  alters  eq.(l7)  t>y  replacing 
p  on  the  right  hand  side  by  p  +  q  as  was  done  in  Ref.  1.  If  one  does 
not  bother  to  center  the  right  hand  side  of  eq.(l5)  in  time  but 
replaces  it  by  its  value  at  t  ,  then  the  transcribed  equation  can  be 
put  into  the  form: 


,n  +  1 


i  + 


1/2 


FU 

Ei  + 


1/2 


<*"  ♦  1/2 


n 

qi  + 
=  0 


1/2 


)<v? 


+  i 

i  +  1/2 


-*”  + 


1/2' 


(18) 


This  equation  does  not  exhibit  the  conservation,  but  it  possesses  it 
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nevertheless  if  one  uses  it  in  conjunction  with  eqs.(lO),  (ll),  and 
(12).  Formally,  the  above  equations  look  like  they  contain  a  vis¬ 
cosity  term  if  one  uses  eq.(l4)  for  q;  and  indeed,  similar  to  the 
viscosity  method  of  von  Neumann  and  Richtmyer  (Ref.  1),  they  lend 
themselves  to  the  treatment  of  not  too  strong  shocks .  On  comparing 
the  results  obtained  by  using  for  q  an  expression  q^  which  is  linear 
in  Au  as  given  by  eq.(l4)  or  an  expression  q^  which  is  quadratic  in 
A  u  as  suggested  in  Ref.  1,  the  following  was  found  (Fig. 4,  p.2l).  The 
use  of  q^  gave  a  fairly  large  overshoot  behind  the  shock  which  damped 
out  quite  rapidly  in  the  wake  of  the  shock.  The  use  of  qg  gave  a 
smaller  initial  overshoot,  but  it  damped  out  slower.  It  seemed 
advisable,  therefore,  to  use  a  combination  of  q^  and  q^  in  order  to 
get  both  a  small  overshoot  and  a  rapid  damping.  Either  one  does 
some  violence  to  the  equations,  which  has  the  effect  of  smearing  out 
discontinuities  such  as  shockfronts,  interfaces,  or  heads  of  rare¬ 
faction  waves  and  of  introducing  disturbances  at  boundaries  which  will 
be  discussed  in  more  detail  later  on  in  this  paper. 

It  is  therefore  desirable  to  keep  this  violence  at  a  mini¬ 
mum,  and  this  may  be  done  in  a  number  of  ways .  The  q ^  suggested  by 
von  Neumann  and  Richtmyer  is  of  the  form 


where 


c  Au  |^u 
2  V 


(19) 


Au 


"  ui  +  1/2 


Ui  -  1/2  * 


(20) 


-11- 


It  was  recognized  quite  long  ago  by  the  people  who  used  this  method 
that  there  was  no  need  for  a  viscosity  term  for  positive  Au,  and  one 


often  finds  it  replaced  by  zero  for  that  range.  This  "spliced  q2" 
has  a  continuous  first  derivative.  One  can  also  adjust  the  parameter 
c2;  and,  in  the  absence  of  q^,  c^  =  2  is  about  right  to  handle  most 
cases.  Of  the  two  expressions  given  for  q^  the  one  in  eq.(l3)  is, 
in  the  absence  of  discontinuities ,  correct  to  one  order  higher  in  &  t 
than  the  one  given  by  eq.(l4).  One  can  reduce  the  violence  done  by 
applying  eq.(l4)  alone  by  mixing  the  two  expressions,  i.e.. 


with  the  indices  as  in  eqs.(l3)  and  (l4).  Aside  from  method  1  , 
which  consists  of  the  continuation  of  q.^  given  by  eq.(2l),  and  the 
spliced  q2,  we  shall  also  test  "method  2",  which  comes  closer  to  the 
requirement  of  having  no  viscosity  for  positive  Au.  In  method  2,  we 
use  q1  +  q2  ,  with  qg  given  by  eq.(l9),  up  to  the  value  of  Au  where 

reaches  its  minimum  and  that  minimum  for  values  of  Au  beyond. 
In  this  method,  q  +  q£  is  also  a  smooth  function  of  Au. 


3 .  Analytical  Solution  of  Test  Problems 


The  above  methods  were  tested  on  a  gas  for  which 


E  =  p  v/U-  1) 


(22) 


with  =  1.4.  The  specific  cases  to  which  tests  were  applied  were  as 
follows : 

A.  A  shock  which  is  induced  by  a  piston  of  constant  velo¬ 
city  uq  running  into  material  of  constant  V  =  ^  and  at  rest.  The 

theory  leads  to  a  shock  with  constant  values  of  u  =  u  „  =  u  , 

V1  do 

V  =  P  =  P2  =  ^  between  the  piston  and  a  shock  front  going 

y+ 1 

with  the  constant  velocity  D.  We  obtain  (with =  -y — j) 


uv 


PlVl 


(>t  -  1) 


(j  -  l)2 


(23) 


and 


"I 


^  +  1 
"  X  + 


D  Z'3' pl^l  •  +  x+  1  ^  ^ 


(24) 

(25) 


A  few  useful  numerical  values  pertaining  to  the  actual  test 
problems  are  collected  in  Table  1. 
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TABLE  1 

Shock  Relations  for  *-=  1.4  (Cases  A  and  B) 


■  ■ 

1 

n 

d/  /p^ 

/7Tn 

4 

21.303 

4.718 

2.427 

5.076 

n.86 

16 

309.4 

5.889 

25.88 

19.27 

50.50 

3.862 

20 

4.654 

2.323 

4.919 

11.41 

16.275 

320 

5.893 

26.70 

19.60 

51.38 
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B.  Same  as  but  with  the  piston  driven  by  a  constant 
pressure.  Tbe  theory  is  the  same  as  under  A. 

C.  A  rarefaction  wave  induced  by  a  piston  receding  with  a 

negative  velocity  uq.  Ibis  leads  to  a  simple  wave  as  discussed, for 

example,  in  Ref.  k,  para.  40.  In  this  solution  one  has  a  fan  in  the 

x  -  t  plane  separating  two  regions  of  uniform  conditions. 

t 


x 

Fig.  2 


(26) 

(27) 


In  the  fan  we  have 


(28) 
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and, with  pV  *  constant  as  in  eq.  (27) > 


*A  -)r+ 1  .  (29) 

/°x  U  Pl  v^y 


The  speed  of  the  front  and  rear  end  of  the  fan  are  obtained  from  (28) 
by  substituting  p.^,  and  po,  Vq  respectively.  Some  numerical  values 
are  found  in  Table  2. 

D.  The  shock  tube  problem  (see  Ref.  4,  para.  8o).  Each  of 
two  media  in  contact  is  initially  at  rest  and  at  constant  density  and 
pressure,  but  the  two  densities  and  pressures  are  different.  Fig.  3 
gives  the  x,t  diagram. 


t 


Between  regions  2  and  3  is  what  Ref.  4  calls  a  contact  discontinuity. 


P2  V1 

itg,,  p  and  u  are  continuous  but  V  is  not.  We  define  ~  f 


and  X  =  Pc/P]/  x  is  8iven»  5  and  are  to  be  found*  We  can  exPress 


U  2  " 


(/<--  1) 


+  1 


u  2 /F  f 
3  =  7-  1  » 


p5  V5 


U-1)  , 

i  -  (|)  2  r 


(30) 

(3D 
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TABLE  2 


Rarefaction  Wave  Relations  for  7"  =  1.4  (Case  C) 


u/  ^P]7! 

■7 

D 

-  1 

.274 

.396 

.390 

-  4 

.000374 

.00356 

.00136 
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and  find  \  "by  equation  u2  =  u^.  ^  and  the  shock  velocity  D  can  then 

be  found  from  eqs.  (2k)  and  (25)  and  v^  from 


v3  =  V^X/S)1^  . 


(32) 


Conditions  in  the  fan  are  given  by  relations 
p  ,  v^,  and  x  replaced  by  p^,  v^,  and  -x. 
are  given  in  Table  3* 


like  (28)  and  (29)  with 
Some  numerical  values 


4.  Numerical  Tests 


The  method  of  integration  was  tested  on  the  four  cases  treated 
analytically  in  the  previous  section.  Equations  (15)  and  (l6)  were 
used  for  conservation  of  mass.  Expressions  (ll)  and  (12)  are  entered 
into  eq.  (10)  to  give  what  was  used  for  the  conservation  of  momentum: 


n+1 

u. 

l 


n 

u.  + 

l 


n  n  n 

g(pi-l/2  *  qi-l/2  ~  Pi+l/2 


PoTxi+l 


qi+l/2^  5t 


(33) 


The  equation  of  state  (22)  is  entered  into  eq.  (l8)  to  give  the 
equations  which  were  used  for  the  conservation  of  energy: 
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TABLE  3 


Shock  Tube  Relations  for  7T  =  1.4  (Case  D) 


V/*  i 

p3/p5 

/3/f  5 

f 2 5 

u3//plVl 

D  llv^l 

2 

.7009 

.7758 

.6357 

.2929 

1.372 

l6 

.2142 

.3328 

.1430 

1.082 

2.077 

32 

.1381 

.2431 

.08251 

1.458 

2.345 

64 

.08711 

.1749 

.04650 

1.743 

2.625 

128 

.05368 

•  1259 

.02563 

2.020 

2.906 

512 

.01923 

.05946 

.007404 

2.552 

3.466 
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where 


n+1 

pi+l/2 


=  y 


n 


Pi+l/2 


+  (  r-  l)(y  -  l) 


/  n 

(pi+l/2 


n 


qi+l/2‘ 


(34) 


y 


vn  / 

i+l/2 

vn+1/ 

i+l/2 


(35) 


Hie  tests  were  carried  out  with  t~  =  1.4.  Runs  were  made  in  the  one- 
medium  cases  A,  B  and  C,  with  30  and  in  the  two-medium  case  D  with 
6o  mass  points. 

In  the  cases  A  and  B,  which  test  the  method  for  its  treat¬ 
ment  of  shocks,  the  problems  were  run  without  interruption  up  to  a 
time  where  one  can  still  run  32  more  cycles  before  the  shock  reaches 
the  last  mass  point.  At  this  time  a  printout  was  made  of  all  physical 
variables  and  in  the  remaining  32  cycles  a  fluctuation  calculation  is 
carried  out.  First,  the  shocked  region  is  located  by  going  from  the 
last  mass  point  backwards  and  testing  for  the  first  pressure  maximum. 

In  the  two  runs  presented  in  Fig.  4,  for  example,  these  are  the  en¬ 
circled  points.  From  there  back  to  the  piston, the  average  pressure  p 
as  well  as  the  maximum  and  the  mean  value  of  the  square  deviation 
A p  =  (p  -  p)  are  calculated  in  each  cycle.  Following  this  through 

—  —  —  o  2  ~2 

32  cycles  the  time  averages  p,  (p  -  p)  ,  Zip  ,  and  Zip  max  were  calculated 

2 

and  the  absolute  maximum  Ap  was  recorded.  Let  us  define 

max  max 

relative  errors  as  follows: 
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Fig.  U  Typical  shock  profiles  with  only  one  type  of  viscosity. 


In  Fig.5>  these  five  relative  errors  are  plotted  as  a  function  of  £ t  for 
a  typical  case.  Also  indicated  is  the  Courant -Friedrichs -Levy  limit, 
which  is  given  by  the  relation  ^  r^p 6  tCFL/  6  x  =  1.  Below 
about  — g -  ,  is  exceedingly  small  and  positive  or  negative  with¬ 

out  any  apparent  correlation  to  any  other  factors.  'Hie  quantities 
^  2  >  ^  3  >an<^  ^ I4.  ^ave  constant  ratios  for  all  cases  considered,  except 
for  small  random  variations.  £,./  £ ^  is  generally  somewhat  larger  for 
small  values  of  3t.  This  is,  however,  not  too  significant;  and  for 
judging  the  merit  of  different  combinations  of  and  c^,we  shall  con¬ 
sider  only  one  of  the  five  E’s ,  namely 

One  can  obtain  some  feeling  for  the  merit  of  different  com¬ 


binations  of  and  in  calculating  this  type  of  flow  by  constructing 
lines  of  constant  in  the  plane.  Such  diagrams  are  presented 

for  two  shocks  of  different  strength  in  Fig.  6.  Optimal  smoothing 


occurs  for  =  1  and  for  values  of  Cg  ranging  about  from  1  to  2.  If 
optimal  smoothing  were  the  only  criterion  for  choosing  the  parameters 


a1  and  the  problem  would  be  solved  with  the  above  observation.  How¬ 
ever,  too  much  smoothing  has  also  some  undesirable  features.  Two  of 
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Fig.  5  Relative  errors  vs.  £t. 
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these  can  be  observed  in  the  above  runs.  Too  much  smoothing  videns  the 
shock  region  and  it  causes  sin  overshooting  of  the  entropy  and  therefore 
too  low  a  density  in  the  zone  next  to  the  piston.  The  origin  of  the 
first  effect  is  obvious,  and  a  numerical  discussion  for  the  quadratic 
viscosity  is  given  in  Ref.  1.  Figure  7  shows  how  this  effect  depends 
on  c^  and  c^  for  two  piston  velocities.  No  special  graph  is  made 
for  the  pressure  boundary  condition  (case  B)  because  it  gives  the  same 
result  as  comparable  cases  A.  The  second  effect  can  be  understood  as 
follows .  After  a  shock  has  travelled  some  distance  away  from  the 
piston,  it  is  essentially  a  steady  state  type  solution  as  discussed  in 
Ref.  1.  One  can  show  that  the  computed  values  of  q  and  v  are  indeed 
close  to  the  theoretical  curve  given  by  eq.  (28)  of  Ref.  1: 

qV  =  UL ±  il  M^(Vi  -  V)(V  -  Vf ).  (36) 

This  relation  is  established  irrespective  of  the  viscosity  law,  pro¬ 
vided  only  that  there  is  enough  viscosity  to  establish  a  steady  state 
solution.  In  the  first  zone,  on  the  other  hand,  one  clearly  has  not 
the  steady  state  type  solution  which  gets  established  after  the  shock 
has  run  some  distance  away  from  the  piston.  In  Fig.  8,  the  dependence 
of  qV  on  V,  as  taken  from  an  actual  run,  is  plotted  both  for  the 
first  zone  and  for  the  interior  zones  and  is  compared  to  the  theoretical 
expression  of  eq.  (36). 

If  we  rewrite  eq.  (18)  in  differential  form 
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Shock  width  vs 


qV/u 


Fig.  8  Comparison  of  viscosity  with  steady  state  theory. 
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dE  +  p  dV 


-  q  dV 


(37) 


and  multiply  by  the  integrating  factor  V**-1,  we  obtain 

4(pv*)  =  -  (2-  i)y  qV^-1  dV.  (38) 

This  relates  q  in  a  simple  fashion  to  the  entropy  increase,  and  an 
inspection  of  Fig.  8  makes  it  apparent  -why  the  entropy  in  the  first  zone 
is  too  large.  It  is  also  apparent  that  the  discrepancy  grows  with  the 
amount  of  viscosity.  Some  measure  of  this  discrepancy  is  the  ratio  of 
the  maximum  of  qV  in  the  first  zone,  which  occurs  at  the  very  start, 
i.e.,  for  V  =  V^,  to  the  theoretical  maximum  for  the  steady  state 
solution  as  calculated  from  eq.  (36).  For  small  St  this  ratio  is 


8 

2  +  1 


u 


i>. 


(39) 


This  of  course  makes  sense  only  for  strong  shocks,  because  the  entropy 
change  for  weak  shocks  is  negligible  anyway.  In  the  case  of  Fig.  8., 
the  contribution  due  to  the  quadratic  term  completely  swamps  the  one 
due  to  the  linear  term.  In  case  B,  i.e.,  when  a  pressure  rather  than  a 
velocity  is  maintained  at  the  piston,  the  depression  of  the  density  of 
the  first  zone  is  less  pronounced  but  nevertheless  present.  Fig.  9  shows 
the  ratio  of  the  density  in  the  first  zone,  as  actually  obtained  in 
some  runs,  to  the  theoretical  density  for  the  steady  state  solution,  as 
function  of  the  parameters  and  c^.  There  are  three  groups  of 
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Fig.  9  Ratio  of  actual  density  in  the  first  zone  to  theoretical 
density  for  the  steady  state  solution. 


-29- 


curves  corresponding  to  three  different  boundary  conditions,  tvo  where 
the  velocity  is  given  and  one  where  the  pressure  is  given. 

Figure  10  shows  some  density  profiles  (skipping  the  first  8  mass 
points)  for  some  selected  combinations  of  c^^  and  Cg.  In  both  the  weaker 
and  the  stronger  shock  curves  the  Courant-Friedrichs-Lewy  number  is  about 
0.4.  A  few  runs  were  also  made  for  infinitely  strong  shock,  i.e.,  shocks 
running  into  a  medium  with  p^  =  0.  The  results  were  in  all  respects  so 
much  like  the  runs  for  which  uq  =  16  y p^  that  it  is  not  necessary 

to  tabulate  or  plot  them.  All  these  data  show  that  for  the  computation 
of  flows  containing' a  shock  one  does  well  using  =  1,  i.e.,  as  large 
as  possible,  and  c^  =  l/2.  This  combination  gives  fairly  good  smoothing 
without  introducing  excessive  troubles  of  the  kind  which  were  discussed 
just  above. 

The  next  point  which  was  studied  was  the  effect  of  viscosity 
pressures  on  the  accuracy  of  calculating  rarefaction  waves  (Case  C). 

Most  calculations  here  were  made  using  method  2  [discussed  just  after 
eq.  (2l)j  .  In  the  calculations  of  shocks  method  1  was  usually  employed 
but  there  is  no  significant  difference  if  one  uses  method  2  instead  be¬ 
cause  the  positive  Au’s  are  not  large  enough  to  reach  the  point  where 

q  is  spliced.  One  can  see  from  eq.  (26)  that  the  fastest  velocity  with 

2  yy  - - 

which  the  material  will  follow  a  receding  piston  is  - — — y  p^  or, 
if  y  =  1.4,  |  UjJ  =  5*916  The  larger  one  chooses  |  uq|  the  more 

severe  a  test  will  one  have  for  the  errors  due  to  the  viscosity  pressure. 
Only  one  run  was  made  with  uq  =  -  ^p^  V^,  for  which  a  pressure 
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DENSITY  PROFILES 


Fig.  10  Density  profiles  for  selected  combinations  of  c^  and  c^. 
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profile  is  presented  in  Fig.  11.  The  other  runs  were  made  with 

uq  =  -  L  ^p^  V-^.  This  is  close  to  the  limiting  case  of  blowing  off  mater¬ 
ial  into  the  vacuum  and  for  that  reason  the  densities  and  pressures  at 
the  "piston"  are  extremely  small.  Fig.  12  presents  a  log  log  plot  of 
the  density  profile  near  the  piston  end  of  the  rarefaction  fan.  Hie 
theoretical  curve  is  obtained  from  eq.  (26)  for  the  plateau  and  from 
eq.  (29)  for  the  fan.  Hie  points  are  obtained  from  a  run  with  q  =  0. 

Hie  differences  here  are  entirely  due  to  the  mesh  size 
=  1,  which  is  large  compared  with  the  Lagrange  distance  x  =  .017 
a  sound  signal  would  have  traveled  on  the  piston  side  of  the  rarefaction 
fan.  If  one  would  plot  the  same  graph  on  a  linear  scale,  the  differ¬ 
ences  would  not  show  up.  The  points  obtained  from  a  run  with  c^  =  1 
and  c  2  =  1  are  so  close  to  the  points  with  c  ^  =  c  2  =  0  actually  drawn 
in  Fig.  12,  that  one  cannot  tell  them  apart.  Hie  difference  between 
runs  with  different  q  does  not  become  significant  until  one  reaches 
points  near  the  head  of  the  rarefaction  wave.  Fig.  13  compares  some 
density  profiles  in  that  range  with  the  theoretical  curve  and  with  the 
curve  computed  with  no  viscosity.  Hie  agreement  is  now  poorer  for 
larger  c,.  It  seems  at  first  surprising  that  the  agreement  gets  better 
for  larger  Cg.  This  is,  however,  easily  understood  if  one  considers 
that  q  is  levelled  off  for  smaller  Au  and  at  a  less  negative  value  for 
larger c g  because  of  the  use  of  method  2  (see  Fig.  Ik). 

Hie  last  group  of  problems  tested  the  shock  tube  problem 
(case  D).  Hie  initial  density  and  pressure  ratio  was  varied  from  2:1 
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Fig.  12  Density  near  piston  end  of  rarefaction  fan. 


Fig.  13  Density  near  head  of  rarefaction  wave. 
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q  BY  METHOD  2 


Fig.  lU  Values  of  q  determined  by  method  2. 
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to  512:1.  The  runs  vith  a  ratio  of  2:1  were  exceedingly  close  to  the 
theory.  Figure  15  shows  a  fairly  typical  density  profile  for  a  ratio 
of  32:1  compared  vith  the  theory.  The  quality  of  the  approximation  in 
the  rarefaction  wave  agrees  with  the  observations  made  on  the  case  C 
runs.  Figure  16  shows  some  runs  for  a  ratio  of  16:1  in  and  near  the 
region  of  the  density  plateau  (region  3  and  parts  of  region  4  of  Fig*  3)* 
The  plateau  shows  up  much  better  in  the  runs  made  with  method  2.  In 
the  shock  region  (region  2  of  Fig.  3)>  however,  method  1  is  better,  as 
evidence  by  Fig.  17,  which  presents  the  error  of  the  average  density 
and  the  average  error  of  the  density  in  that  region. 


-37- 


-38- 


Fig.  16  Density  near  plateau. 
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Fig.  17  Error  of  the  average  density  and  the  average  error  of  the 
density  in  the  shock  region. 
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